Table 1 presents the details of RECL and OWA values (given to the nearest one-tenth millimeter) used in the present study for each fossil hominin specimen, as well as the mean values for fossil and extant species.

Most of the micro-ct measurements used in this study were employed previously by Braga et al. (2015) (see Table S1 for more details) with isometric voxel dimensions ranging from 7.0 to 41 microns (μm). Therefore, we do not repeat here the details about the samples and the body mass data already given in Braga et al. (2015). Additional measurements of RECL and OWA were taken from the micro-cts of 7 Gorilla gorilla skulls (from specimens housed at the “Musée royal de l’Afrique centrale”, Tervuren, Belgium; with only one of them of known sex) obtained using the Nanotom (GE Sensing) at the “Fédération de recherche” FERMAT (CIRIMAT, Toulouse) (all with isometric voxel dimensions of 32.6 μm). The micro-ct data for StW 498e, StW 151 and StW 53 were obtained in South Africa using the XTH 225/320 LC dual source system (Nikon) at the Palaeosciences Centre, University of the Witwatersrand, Johannesburg, with isometric voxel dimensions ranging from 28.1 to 30.5 μm.

The micro-ct measurement methodology of RECL and OWA is detailed in Braga et al. (2015). Moreover, we illustrate the micro-ct measurement protocol used in this study for OWA with an example taken for the P. robustus SK 879 specimen (Fig. 2). This example is motivated in part by the difficulties expressed in the literature about accurate measures of the oval window in the narrow oval window niche (fenestra vertibuli), a variably deep depression located partly behind the overhanging and prominent promontory of the tympanic cavity (promontorium tympani) (Fig. 2). In a detailed analysis of the stapes footplate, Sim et al. (2013) stated that measurements taken on photographs “might not be the most precise way of measuring the dimensions of the ossicles” (op. cit., p. 160). Noussios et al. (2016) also considered that photographic measurements of the stapes footplate in the oval window niche may entail errors regarding the projection of such a small surface into a screen (parallax error). Since several studies now consider in situ micro-ct measurements of the middle ear as more accurate and efficient than measurements taken on photographs (this point will be discussed in more details below), we prioritize micro-cts to measure OWA in fossil hominins and extant catarrhines (Braga et al., 2013 and Braga et al., 2015).

We first visualized (by using the Avizo solftware. https://www.fei.com/software/amira-avizo/) the oval window niche in 3D after extracting an isosurface from the micro-ct data set, as shown in Fig. 2a with the case of SK 879. We then defined an oblique slice that best-fitted the complete outline of the oval window at the location of the annular stapedial ligament, a ring of fibrous tissue that connects the oval window to the stapes footplate. The OWA was then measured from its segmentation on this oblique slice. Since we used high resolution micro-ct data (Table 1), the outline of the insertion area of the annular stapedial ligament could be visualized accurately. This allowed us to reduce to a minimum (less than 5%) the interobserver variability between two operators (J.B. and J.R.D.) in setting an oblique slice in the plane of the OWA.

All the phylogenetic methods and tests presented below have been implemented by using the statistical maximum likelihood (ML) framework that was first introduced for phylogenetic studies by Edwards and Cavalli-Sforza (1964) for gene frequency data. The application of ML to phylogenies involves searching for the single tree (i.e. a topology with branch lengths) that, under a given evolutionary model, best explains the observed data (i.e. the traits for each species, for more details, see Baum and Smith, 2013). In doing so, it is conventional to record the logarithm of the likelihood, the log-likelihood (to avoid problems associated with handling very small numbers) of that tree as the main criterion.

An important point to stress here is that analyses for morphological data are still often based on parsimony, a method that is less computationally demanding and “consistent if the rate of evolutionary change per branch of the tree is sufficiently small” (Felsenstein, 2004:122). Intuitively speaking, the statistical property of parsimony (for more details, see Felsenstein, 2004) means that it minimizes the number of changes that have occurred in each branch and therefore assumes a low rate of change (i.e. fewer evolutionary events, so the phylogram will show relatively shorter branch lengths in all characters). Moreover, the advantage of ML methods over parsimony analysis is that the former use information on branches lengths available in phylograms. Given the fact that different characters are often expected to evolve at different rates and that longer branches offer more opportunities for changes to occur than shorter branches, ML methods are particularly useful in this context. We use freely available packages (listed below) in the R Statistical Computing environment http://www.r-project.org/. These packages are summarized in the CRAN project: http://www.cran.r-project.org/web/views/Phylogenetics.html.

We first investigate whether the cochlear data observed for extant species derive from a nonrandom, gene-based tree-like descent process. The basic principle is to determine whether a given tree derived from genetic distances better fits a set of cochlear data (as represented by species mean values) at its tips, when compared with the fit obtained when the cochlear data have been randomly permuted across the tips of the tree, thus destroying any phylogenetic signal that may have existed. We investigate the phylogenetic signal in this dataset before and after controlling for body size. We use distinct R packages to compare two estimates of the phylogenetic signal for RECL and OWA considered separately: the λ of Pagel (1999) the K of Blomberg (Blomberg et al., 2003).

Pagel's (1999) λ is a parameter multiplied to each of the off-diagonal elements of the variance-covariance matrix (see Nunn, 2011, for more details). When λ = 0, the tree has a single polytomy at the basal node for all species, whereas when λ = 1 the original candidate tree is recovered. Statistical tests for phylogenetic signal are performed under the null hypothesis that λ = 0. Tests for less signal than the candidate tree are performed under the null hypothesis that λ = 1. We compute Pagel's λ with the R package “Geiger” (its function “fitContinuous”) under the assumption of a Brownian model of evolution (see below). We also compute Pagel's λ with the R package “Caper” (its function “pgls”) with no prior assumption about any given model of evolution. In all instances, we test our null hypotheses with the likelihood ratio (LR) test (twice the log-likelihood difference between the two models is expected to fit a Chi² distribution with P degrees of freedom) (Baum and Smith, 2013 and Nunn, 2011).

We compute Blomberg's K statistic (Blomberg et al., 2003) with the R “Picante” and “Ape” packages in order to test also whether the observed distributions for RECL and OWA across catarrhine species exhibit more or less divergence than expected when evolving under Brownian motion. Values of K range from 0 to infinity, with K = 1 indicating Brownian motion evolution. K > 1 indicates that close relatives are more similar than expected, and K < 1 indicates more divergence between taxa than expected under a Brownian model. In order to decide whether the calculated K value corresponds to a Brownian model, we test if it is significantly higher than the mean K value obtained after simulating 1000 datasets (i.e. trees) with random permutations of the tips.

We also use phylogenetic generalized least squares (PGLS) (with the “pgls” function of the R package “Caper”) to assess whether either RECL or OWA considered separately may have undergone correlated evolution with body mass (BM). Therefore, we estimate the slopes of RECL versus BM and OWA versus BM regressions, and test whether these slopes are significantly different from 0. We subsequently calculate the phylogenetic signal of RECL and OWA after allometric correction (i.e. removing the allometric effect of body mass here represented by mean species values).

In this part of the study, we aimed to characterize the tempo (or rate) and the mode (the mechanisms) of RECL and OWA evolutionary changes by testing which evolutionary model best explains the observed data. We consider the following three evolutionary models: BM, OU and ACDC.

We first consider the Brownian motion (BM) model assuming that rates of evolution are constant over time (i.e. larger changes more likely occur on longer branches). Therefore, this model was also called “constant variance process” (Freckleton et al., 2002 and Pagel, 2002). In this model, the instantaneous phenotypic variance σ ² per unit time t (i.e. the change along a branch of length t, or rate of change) is drawn from a normal distribution with mean 0 and variance σ ² t. Therefore, the trait variance is proportional to time (or branch length) and the displacements in different branches of a tree are independent. It is important to add here that the BM model does not necessarily indicate neutral or random evolution. It can instead reflect adaptive evolutionary changes (Nunn, 2011).

In contrast to BM, the Ornstein–Uhlenbeck (OU) model can be viewed as acting to limit evolutionary variation. Indeed, it is biologically unrealistic to consider that the range of a trait value could vary infinitely. The OU model can therefore be applicable when traits likely evolved with constraints on their maximum or minimum values, depending on the strength of stabilizing selection (Felsenstein, 1988 and Garland et al., 1993). It is a BM model pulled to include one or more selective optima that assert an attractive force on Brownian trait evolution. When the strength of this attraction is zero, the OU model is identical to a BM model. Very strong stabilizing selection can obliterate the effects of history such that phylogenetic signal disappears. In this event, trait values will be more similar in closely or distantly related species experiencing the same selective pressures.

The acceleration or deceleration (ACDC) of Brownian motion evolution from the root to the tips of the tree is the third model used in this study. The ACDC model (also called “early burst”) describes evolution that either increases (accelerates, AC) or decreases (decelerates, DC) in rate over time.

We choose the evolutionary model that best describes the species data sampled in our study with the use of LR tests and the Akaike information criterion (AIC) (Akaike, 1974). The model with the smallest AIC is the preferred one. We also use the corrected AIC (AICc) for small samples.

We estimate RECL and OWA ancestral values (using the preferred evolutionary model chosen in the previous step) at internal nodes of a phylogeny reconstructed with morphological and molecular data by grafting a tree for three extinct hominin species of known cochlear values (A. africanus, P. robustus and H. erectus) onto the gene-based phylogram representing contemporaneous catarrhine species (Fig. 1). The extinct hominin phylogram (i.e. branch lengths proportional to the amount of evolutionary change) is obtained from Organ et al. (2011) and was inferred using 109 morphological characters from Strait and Grine (2004) in a Bayesian framework.

We use the R package “Phytools” (and its “fastAnc” function) (Revell, 2012) to compute ancestral values and 95% CIs (to reflect uncertainty) using ML (i.e. finding values that maximize the probability of our data). We consider that our approach is unbiased because we use the evolutionary model that best fits our data. The ability of different methods (e.g., ML under different models of evolution) to accurately estimate ancestral values and their associated errors has been investigated by Martins (1999), who showed that a BM model performed reasonably well, especially for more recent ancestors. Therefore, in order to reduce uncertainty, we limit our estimates for hominin ancestral nodes that are close to the tips (as opposed to those that are deeper in the tree). We consider the following four nodes: (Pan, Hominins), (A. africanus, P. robustus, Homo), (P. robustus, Homo) and (H. erectus, H. sapiens) (Fig. 1).

We first represent in a simple 2D-point graph, the variation of RECL versus OWA among early hominins, H. sapiens sapiens, Gorilla gorilla and Pan troglodytes/paniscus (Fig. 3). At first glance, we observe a clear separation between a group comprising the four Australopithecus specimens from Sterkfontein (StW 498e, StW 255/259, StW 98 and StW 329), StW 53, StW 151 and KB 6067, and a second group made up of P. robustus from Swartkrans (SK 879 and SkW 18) and H. erectus (SK 847) with noticeably higher values for both RECL and OWA (Fig. 3). In this regard, it is important to note that RECL is a better discriminator than OWA, the SK 847 OWA value falling within the variability sampled for specimens derived from Sterkfontein Member 4. We also observe that H. sapiens and Gorilla gorilla specimens overlap one another and that most of them have higher values of RECL than that for Pan troglodytes/paniscus. While the three specimens of P. robustus from Swartkrans and H. erectus group fall within H. sapiens and Gorilla gorilla ranges, all the other fossil hominins sampled here from Kromdraai and Sterkfontein group with values for Pan troglodytes/paniscus. Moreover, the StW 329 specimen shows remarkably low RECL and OWA values.

For both RECL and OWA, we obtain a high Pagel's λ that is significantly different from 0 but not different from 1, either under the assumption of a Brownian model of evolution, or with no prior assumption about any given model of evolution (Table 2). The same trend is observed when we compute the Blomberg's K statistic. In all instances, the K values are significantly higher than those obtained randomly and indicate either Brownian motion evolution (RECL; K = 0.98) or closer relatives more similar than expected (OWA; K = 1.52) (Table 2).

Since BM shows a high phylogenetic signal for both Pagel's λ (significantly different from 0, but not from 1) and Blomberg's K (1.22) (Table 2), it is important to investigate whether RECL and OWA still show a high phylogenetic signal after controlling for BM. We obtain a low Pagel's λ (0.32; not significantly different from 0) for RECL after controlling for BM (r ² = 0.78, P < 0.005). For OWA, after controlling for BM (r ² = 0.68; P < 0.005), the Pagel's λ is higher than for RECL, with a value (0.735) significantly different from 0, but not different from 1 at 0.03% (Table 2).

We first compute the log-likelihood for a Brownian model of evolution for RECL and OWA considered separately (Table 3). We then compute the log-likelihoods for OU and ACDC models of evolution of RECL and OWA in order to determine whether one of them better explains our observed data than Brownian motion. Both LR tests and AIC criteria (either AIC or AICc) show that the evolutionary model that best explains the data (i.e. RECL and OWA values and the gene-based phylogram) observed in extant species is Brownian motion (Table 3). This model is therefore used for ancestral states reconstructions.

Since both RECL and OWA show high phylogenetic signals (when considered separately), we estimate the ancestral values (with 95% CIs) of these two parameters at the three hominin nodes in the phylogram (Fig. 1), as well as at the Pan-Hominin node (Fig. 4 and Fig. 5).

With regard to RECL, there is a clear trend for an increase of the ancestral values from the oldest node (Pan, Hominins) to the youngest one and (H. erectus, H. sapiens) (Fig. 4). The (A. africanus, P. robustus, Homo) node's CI at 95% fits completely within the (Pan, Hominins) node's CI at 95%. Moreover, the (P. robustus, Homo) and (H. erectus, H. sapiens) nodes’ RECL estimations fall outside the (A. africanus, P. robustus, Homo) node's CI at 63% (i.e. with only one standard deviation) and 95%, respectively. There is no overlap between the (A. africanus, P. robustus, Homo) and the (H. erectus, H. sapiens) nodes’ CIs at 63% (Fig. 4).

There is also a trend for an increase of the OWA ancestral values from the oldest node (Pan, Hominins) to the youngest one (H. erectus, H. sapiens) (Fig. 5). As for RECL, the (A. africanus, P. robustus, Homo) node's CI at 95% is encompassed by the (Pan, Hominins) node's CI at 95%. However, there is more overlap between the three OWA hominin nodes’ CIs at 63% (Fig. 5) than for RECL (Fig. 4). In particular, there is an almost complete overlap between the (A. africanus, P. robustus, Homo) and the (H. erectus, H. sapiens) nodes’ CIs at 63% (Fig. 5).

In order to evaluate the impact on estimates of RECL and OWA ancestral values when grafting a tree for three extinct hominin species onto the gene-based phylogram representing extant species only, we computed the (Pan, Homo) node CIs before the incorporation of fossil data (as in Braga et al., 2015). For both RECL and OWA, the (Pan, Homo) nodes’ CIs at 95% (12.5–14.7 mm and 2.6–3.9 mm², respectively) were almost identical to the (H. erectus, H. sapiens) nodes’ CIs obtained from the tree combining extant and fossil species (Fig. 4 and Fig. 5). This result demonstrates the usefulness of grafting fossil tip data onto the gene-based phylogenetic tree for ancestral reconstructions at some nodes.

Apart from Braga et al., 2013 and Braga et al., 2015, we could not find any other RECL measurements in early hominins to compare with ours. Moreover, we discuss below why we consider that the micro-ct OWA measurements presented here are not interchangeable with published measurements of the oval window or stapes footplate area taken on scaled digital images in fossil hominins (e.g., Quam et al., 2015). Quam et al. (2015) assessed photographically the OWA or the stapes footplate area and considered that, in all instances (i.e. for all fossil specimens), the latter represented 90% of the former. If we use this 90% proportion to compare the micro-ct OWA values in fossil hominins presented in this study (Table 1) and the photographic measurements published in Quam et al. (2015, Table 1) (e.g., the stapes footplate area of SKW 18 or the OWA of SK 879), we observe marked discrepancies (except for STW 328), in particular for P. robustus. Here we take one example presented in Fig. 2 in order to briefly discuss this point. The most important discrepancy (1.6 mm²) is observed for the SK 879 specimen illustrated with an oblique slice that best-fits the complete outline of its OWA (Fig. 2B) and with a 3D view of its oval window niche (Fig. 2C). As shown in a comparison between Fig. 2B and C, the configuration of the oval window niche poses two main difficulties when trying to measure the small OWA with no parallax error.

The first difficulty is due to the fact that, even when the middle ear cavity is visually accessible without endoscopy (as for SK 879), the promontory of the tympanic cavity prevents the visualization of the complete outline of the oval window to allow for an accurate photographic measurement of OWA to the nearest one-hundredth square millimeter with no potential parallax error. This potential parallax error explains why in situ micro-ct measurements of middle ear features have been considered as more accurate than photographic measurements (e.g., Noussios et al., 2016 and Sim et al., 2013).

The second difficulty is related to the need for a perfect calibration of the photograph to be used for accurate measurements of OWA. First, the scale used for calibration needs to be placed at precisely the same distance from the camera lens as the oval window itself. Moreover, if a millimeter ruler scale is used, it needs to fit inside the narrow oval window niche next to the OWA (Fig. 2C), and also to lie exactly in its plane (i.e. the OWA and the ruler scale need to be co-planar). Small misalignment and/or placement errors of the ruler scale will result in uncertainties on the calibration of the photograph used for measurements.

In the case of the SK 879 micro-ct data, we obtained a 4.3-mm² OWA value (Table 1). In Fig. 2B, we superimpose a grid of 0.5 × 0.5 mm uniformly distributed squares (with surface areas of 0.25 mm² each) over the oblique slice that we used to obtain this OWA value for SK 879. Fig. 2C shows exactly the same viewpoint, magnification and orientation of this specimen as in Fig. 2B. The same grid is superimposed on a 3D endo-tympanal view of the SK 879 oval window niche. A visual comparison shows a noticeable difference between Fig. 2B and C in the measure of the OWA with at least 0.5 mm² (2 squares). As shown here, such a difference is due to the overhanging promontory of the tympanic cavity that could be also seen in a photograph or under a microscope. Moreover, the grid superimposed on the OWA of the SK 879 specimen (Fig. 2B) allows the reader to visually assess that its value is close to 4.0 mm². A 2.43-mm² value was given in Quam et al. (2015) for the SK 879 estimated stapes footplate area. If we assume that the stapes footplate represented 90% of the OWA in SK 879 (as suggested in Quam et al., 2015), we obtain a 2.7-mm² value. This measurement is incompatible with the micro-ct value (4.3 mm²) that we repeatedly obtained for SK 879 after careful segmentations (Fig. 2).

Another important point to discuss is the correspondence between the stapes footplate area and the OWA. The difference between these two areas is due to the size of the annular ligament of the stapes. In a detailed micro-CT measurement study of the annular ligament of the human stapes, Mohammadi et al. (2016) found it highly variable in its volume and thickness. Because of this variability among modern humans, we would urge caution when estimating OWA from the stapes footplate area in other extant species, as well as in fossil hominins with a single constant proportion. The proportion of 90% used for fossil hominins in Quam et al. (2015) likely represents only the upper part of the variation of the relationship between OWA and the stapes footplate area. Here, we argue that more comparative micro-CT measurements of the OWA and the stapes footplate area are needed in both modern humans and other extant catarrhines for better estimations of measurements errors when assessing the former variable from the latter. Indeed, from the micro-ct data of the middle ear obtained by one of us (J.B.), we observe a 76–87% variability (caused by interindividual differences and lateral asymmetry) in the proportion of the OWA represented by stapes footplate area in only 12 measurements for modern humans. This observation is well in line with the results presented in Mohammadi et al. (2016). Moreover, when we investigate the proportion of the OWA represented by the stapes footplate area using our micro-ct data for 11 specimens representing non-human primate species, we observe a range between 79% (Papio hamaryas) and 90% (Pongo pygmaeus).

The present study confirms, using additional tests, the very high phylogenetic signal of RECL and OWA parameters and their close correlation with body mass, as already observed by Braga et al. (2015). Moreover, we newly observe that a Brownian motion model is the best evolutionary framework to compare RECL and OWA differences between catarrhine species, including fossil hominins. Importantly, while RECL evolution is associated with evolutionary changes in body mass, OWA likely evolved more independently. Therefore, comparisons between fossil hominin specimens and species based on RECL and OWA may lead to reliable identifications of monophyletic groups with their associated synapomorphies (as opposed to homoplasies).

The fossil cochlear values presented here confirm the RECL distinctiveness of the Paranthropus and H. erectus specimens found at Swartkrans compared to the Australopithecus specimens from Sterkfontein (Braga et al., 2013 and Braga et al., 2015). The Paranthropus and H. erectus RECL values are similar to measurements for extant humans and gorillas, but also noticeably higher to those for any specimen from Sterkfontein, including StW 151 or StW 53 (Fig. 3). Interestingly, the StW 53 RECL and OWA values appear very close to those recently published for KB 6067, a juvenile fossil hominin from the South African site of Kromdraai (Braga et al., 2013 and Braga et al., 2016) which has been attributed to an evolutionary less derived P. robustus form than that represented at Swartkrans. Indeed, as noted by Braga et al. (2013), StW 53 and KB 6067 also appear similar in the morphology of their semi-circular canals. Interestingly, while StW 53 and StW 151 petrous bones have been regarded as very close in morphology and more like the modern human condition than exhibited by other Sterkfontein specimens, the OWA size in StW 53 falls in the range of P. paniscus whereas the much larger OWA of StW 151 falls in the range of both H. sapiens and G. gorilla. Moreover, the RECL values of both StW 53 and StW 151 fall well below the P. robustus. H. erectus and H. sapiens values (Fig. 3).

As indicated in Braga et al. (2013), the taxonomic attribution of KB 6067 from Kromdraai Member 3 is as yet unclear. A firmer taxonomic attribution of the oldest Kromdraai paranthropines (including KB 6067) is still pending on the analysis of the larger hominin sample discovered in the oldest fossiliferous deposits from this site (Braga and Thackeray, 2016 and Braga et al., 2017). This will allow us to determine whether the southern African Paranthropus hypodigm represents either a single and variable P. robustus species with a time span yet unknown, or distinct taxa (i.e. at most the more plesiomorphic P. robustus from Kromdraai and P. crassidens from Swartkrans; at least a single evolutionary P. robustus species consisting in a lineage of ancestral descendant populations; for more details, see Braga et al., 2016 and Braga et al., 2017). Further discoveries combined with progress in measuring the complex cochlear shape will allow us to test whether the variable OWA and RECL signatures obtained between Swartkrans (e.g., SKW 18 and SK 879) and Kromdraai (e.g., KB 6067) paranthropines are consistent with “a plesiomorphic status of the KB hominins, indicative of their ancestral status for the P. robustus + boisei clade, giving rise before 2.3 Ma to the split of P. boisei in East Africa and P. robustus survivors in South Africa” (Braga et al., 2016:65).

Even though our results are preliminary, given the phylogenetic tests presented here and their advantage over parsimony-based interpretations (see below), we argue that StW 53 and KB 6067 may represent one or two distinct smaller-bodied, less derived hominin form(s), as compared to Paranthropus and H. erectus specimens represented at Swartkrans. We acknowledge that labyrinthine similarities are not sufficient to establish firm taxonomic affiliations and that our conclusions need to be tested by further comparative morphological studies. For instance, although the spiraled shape of the cochlea appears globally similar among catarrhines, the development of new 3D morphometric methods may uncover taxonomic distinctions because its spiraled structure makes it inefficient to study with standard Euclidian metrics.

The base of the cranium has often been considered evolutionarily more conservative that the face and teeth (Bosma, 1976 and de Beer, 1937; Lieberman et al., 1996, Strait, 1998 and Strait et al., 1997). Variations in the cranial base other than those investigated here (e.g., the angle of the petrous part of the temporal bone to the coronal plane) have been considered useful in phylogenetic analyses (Dean, 1986 and Dean and Wood, 1981), but have not been tested using the statistical methods employed here. In line with our results, these basicranial features have been regarded as distinct between the contemporaneous great apes and Australopithecus, on the one hand, and Paranthropus and Homo on the other. Further comparative data and investigations at both intra- and inter-specific levels in extant and fossil taxa will improve our phylogenetic interpretations, and in particular, we need data on the cochlear variability in highly dimorphic hominid species such as G. gorilla. As yet, we only investigated the RECL variability in a G. gorilla sample of very limited size (n = 14) and with most specimens (n = 12) of unknown sex (n = 12) and two males.

On the basis of the high phylogenetic signal of RECL and OWA and our evolutionary model-based phylogenetic statistics, here we interpret the P. robustus and Homo shared and derived cochlear features reported on Swartkrans specimens (i.e. higher values for both RECL and OWA) as synapomorphies. From our very limited sample of fossil hominins and the results presented in Fig. 4 and Fig. 5, it is tempting to suggest a A. africanus to P. robustus ancestral-descendant relationship. However, such a conclusion would be unwise if it rests only on high RECL and OWA phylogenetic signals. We nevertheless argue that the improvement of our understanding of the phylogenetic relationships between the southern African paranthropines and other early hominin taxa, including early Homo and Paranthropus from East Africa, depends on further analyses of the skull base, including the morphology of the cochlea.

As yet, the cranial base flexion has garnered a lot of attention and, to the best of our knowledge, its derived morphology shared between Homo and Paranthropus has not been explored among fossil and extant primates using phylogenetically informed statistical analyses. Even though we used such analyses in the present cochlear study, we cannot yet determine whether the RECL and OWA derived features shared between P. robustus and Homo at Swartkrans are also shared with eastern African Paranthropus and Homo. There is no consensus on the basal species leading to Paranthropus and Homo, depending on the phylogenetic role accorded to A. africanus and P. aethiopicus (Strait and Grine, 2004) and Paranthropus phylogeny (Wood and Constantino, 2007, fig. 10). Therefore, more comparative studies on Australopithecus and Paranthropus specimens from southern and eastern Africa will be needed to decide whether A. africanus represents the most likely common ancestor of P. robustus and Homo, Paranthropus and Homo or if an alternative scenario is more likely.

In line with results from Baker et al. (2016) and the positive selection associated with shifts in molar area relative to body size in Paranthropus, our phylogenetic findings support a significant role of selection in shaping the evolution of P. robustus and Homo auditory capacities, at least as represented by RECL (Fig. 4). The positive association between BM and RECL observed in this study (see also Braga et al., 2015) means that there is an evolutionary correlation between BM and cochlear size among catarrhines, including fossil hominins. Since RECL does not show a phylogenetic signal significantly different from 0 after controlling for BM, the ancestral values proposed here highlight, at least partly, the role of selection on a larger body size in shaping cochlear evolution in Homo. Only further measurements of RECL in fossil Homo specimens will allow us to determine the exact shifts in cochlear size relative to body size during the earliest part of the evolution of our genus. The use of a phylogram combining extant and fossil species to estimate ancestral RECL values at three hominin nodes in this study (Fig. 4), confirms a previous finding that RECL is “hypertrophied” in the genus Homo, as opposed to the australopiths (Australopithecus and Paranthropus) (Braga et al., 2015). Indeed, when RECL is represented in a phylo-morphospace (Fig. 4), the (H. erectus, H. sapiens) nodes’ value is clearly shifted towards a higher RECL, as compared to the ancestral value at the (P. robustus, Homo) node. Moreover, when we consider the OWA parameter with its high phylogenetic signal after correcting for BM (contrary to RECL), we observe that OWA values are shifted towards higher values at the (P. robustus, Homo) node. Further estimates of ancestral node values will be needed to help refining the predictions by adding new fossil hominin species values at the tips of the tree used in the statistical analyses. To this end, cochlear investigations on A. afarensis, P. aethiopicus and P. boisei will be crucial.

The present study demonstrates that a simple Brownian model of evolution best explains the observed RECL and OWA data in contemporaneous species. Therefore, this model assuming constant rates of evolution over time is appropriate for ancestral reconstructions. Moreover, this result highlights the importance of investigating cochlear changes among fossil hominins by using an evolutionary framework with stretched and compressed branch lengths that conforms to our assumed underlying Brownian motion mode. Such an evolutionary framework is more realistic because it assumes that changes in character values are more likely along longer branches, i.e. over longer stretches of time. When branch lengths are available, our proposed model-based method has the advantage of incorporating these data. By contrast, the parsimony-based methods used for cladistics reconstructions assume that changes are equally likely on all branches of the tree, regardless of their length. This assumption would be safe only for very slowly evolving traits, but not for cochlear evolution among hominins. In this context, methods using parameters with a measured phylogenetic signal (i.e. statistically tested using a phylogenetic comparative approach) will be very useful to investigate further the hypothesized Paranthropus monophyly (see reviews in Wood and Constantino, 2007 and Wood and Boyle, 2016), its diversification into eastern and southern African forms from a common ancestor often set during the 2.7–2.3 Ma period, and its phylogenetic relationships with the origin of our own genus Homo.